Abstract
This notebook is a template for evaluating COVID-19 hospitalization forecast submissions from COVIDhub. After inputting a set of parameters (forecasters, COVID signals, etc), the template yields a comprehensive report on the predictions of COVID forecasters as well as their performance compared to the ground truth. The visualizations generated by the template offer an intuitive way to compare the accuracy of forecasters across all US states.
\[\\[.4in]\]
Every week, forecasters submit their hospitalization predictions to COVID-19 ForecastHub. In this report, we rely on an AWS bucket that contains the estimates of a handful of signals (e.g., COVID death, cases, hospitalization, etc). Furthermore, the AWS server stores an array of evaluation metrics of these forecasts (e.g., Absolute Error, Weighted Interval Score, and 80% Coverage). Alternatively, the data can be retrieved from the publicly accessible covidcast and covideval APIs.
s3a <- get_bucket("forecast-eval", region = "us-east-2")
s3b <- get_bucket("forecasting-team-data")
scores1 <- s3readRDS("score_cards_state_hospitalizations.rds", s3a,
region = "us-east-2")
scores1 <- subset(
scores1,
forecaster %in% union(params$forecasters, "COVIDhub-baseline"))
# The Hub always makes forecasts on Monday (dofw 2)
# We move all forecast_dates to the following Monday and shorten the ahead
wday_shift <- function(x, base_dofw = 2) (base_dofw - wday(x)) %% 7
scores1 <- scores1 %>%
mutate(ahead = ahead - wday_shift(forecast_date),
forecast_date = forecast_date + wday_shift(forecast_date))
# need to pass in the right forecaster here
scores <- s3readRDS(params$aws_scores, s3b) %>%
magrittr::extract2("reformatted.buggy.matched.scorecards") %>%
select(ahead, geo_value, forecaster, forecast_date, target_end_date,
actual, wis, ae, cov_80, value_50, data_source, signal,
incidence_period)
# Logan's names are brutally long...
our_pred_dates <- unique(scores$forecast_date)
forecast_dates <- our_pred_dates
aheads <- unique(scores$ahead)
geo_values <- unique(scores$geo_value)
veggies <- c("asparagus", "broccoli", "chard", "daikon",
"escarole", "fennel", "garlic", "horseradish",
"jicama", "kohlrabi", "lettuce", "mushroom",
"nori", "okra", "parsnip", "radish", "squash",
"tomato", "watercress", "baseline1", "yuca")
names_table <- tibble(forecaster = unique(scores$forecaster), veggies = veggies)
scores <- left_join(scores, names_table) %>%
select(-forecaster) %>%
rename(forecaster = veggies)
results <- scores1 %>%
select(ahead, geo_value, forecaster, forecast_date, target_end_date,
actual, wis, ae, cov_80, value_50, data_source, signal,
incidence_period) %>%
bind_rows(scores) %>%
filter(forecast_date %in% forecast_dates,
ahead %in% aheads,
geo_value %in% geo_values)
The target forecast dates are:
2020-11-16, 2020-11-23, 2020-11-30, 2020-12-07, 2020-12-14, 2020-12-21, 2020-12-28, 2021-01-04, 2021-01-11, 2021-01-18, 2021-01-25, 2021-02-01, 2021-02-15, 2021-02-22, 2021-03-01, 2021-03-08, 2021-03-15, 2021-03-22, 2021-03-29, 2021-04-05, 2021-04-12, 2021-04-19, 2021-05-03, 2021-05-10, 2021-05-17, 2021-05-24, 2021-05-31, 2021-06-07, 2021-06-14, 2021-06-21, 2021-06-28, 2021-07-05, 2021-07-12, 2021-07-19, 2021-07-26, 2021-08-02, 2021-08-09, 2021-08-16, 2021-08-23, 2021-08-30, 2021-09-06, 2021-09-13, 2021-09-20, 2021-09-27
The template will compile data of the following forecasters:
COVIDhub-ensemble, USC-SI_kJalpha, MOBS-GLEAM_COVID, Karlen-pypm, JHUAPL-SLPHospEns.
For this analysis, all of Logan’s forecasters have been renamed:
kableExtra::kbl(names_table) %>%
kableExtra::kable_styling(bootstrap_options = c("striped", "hover"))
| forecaster | veggies |
|---|---|
| ______________________none dacouno zacouno neimacnah | asparagus |
| ______________________none dacouno zacouno seimacnah | broccoli |
| ______________________none oacouno zacouno neimacnah | chard |
| ______________________none oacouno zacouno seimacnah | daikon |
| ______________________none oareigo zareifo neimacnah | escarole |
| ______________________none oareigo zareino neimacnah | fennel |
| ______________________none oareino zareino neimacnah | garlic |
| ______________________none zacougo dacouno neimacnah | horseradish |
| ______________________none zacougo dacouno neiwinnah | jicama |
| ______________________none zacouno dacouno neimacnah | kohlrabi |
| ______________________none zacouno dacouno neiwinnah | lettuce |
| __________________never-up oareigo zareino neimacnah | mushroom |
| _________________always-up oareigo zareino neimacnah | nori |
| ________________always-low oareigo zareino neimacnah | okra |
| _______________always-down oareigo zareino neimacnah | parsnip |
| _____________always-steady oareigo zareino neimacnah | radish |
| ___marginal_cheating_u/sdl oareigo zareino neimacnah | squash |
| __nqcomb_logistic1_u/s/d/l oareigo zareino neimacnah | tomato |
| _marginal_cheating_u/s/d/l oareigo zareino neimacnah | watercress |
| baseline1 | baseline1 |
| marginal_logistic1_u/s/d/l oareigo zareino neimacnah | yuca |
\[\\[.07in]\]
Mean <- function(x) mean(x, na.rm = TRUE)
GeoMean <- function(x, offset = 0) exp(Mean(log(x + offset)))
facets.label = str_glue("{aheads} days ahead")
names(facets.label) = aheads
subtitle = sprintf("Forecasts made over %s to %s",
format(min(forecast_dates), "%B %d, %Y"),
format(max(forecast_dates), "%B %d, %Y"))
plot_wis <-
plot_canonical(results,
x = "forecast_date",
y = "wis",
aggr = GeoMean,
grp_vars = c("forecaster","ahead"),
facet_rows = "ahead", dots = FALSE,
base_forecaster = "COVIDhub-baseline") +
labs(title = subtitle,
x = "Forecast Dates",
y = "Geometric Mean WIS") +
#geom_point(aes(forecast_date, round(wis, digits = 2), color)), alpha = 0.05) +
facet_wrap(~ahead, nrow = 4, labeller = labeller(ahead=facets.label)) +
theme_bw() +
theme(plot.title = element_text(hjust = "center"),
legend.position = "bottom",
legend.title = element_blank()) +
scale_y_log10() +
geom_hline(yintercept = 1, size = 1.5) +
scale_color_viridis_d() +
guides(color = guide_legend(ncol = 2))
ggplotly(plot_wis, tooltip="text", height=800, width= 1000) %>%
layout(hoverlabel = list(bgcolor = "white"))
plot_wis <-
plot_canonical(results,
x = "forecast_date",
y = "wis",
aggr = Mean,
grp_vars = c("forecaster","ahead"),
facet_rows = "ahead", dots = FALSE,
base_forecaster = "COVIDhub-baseline") +
labs(title = subtitle,
x = "Forecast Dates",
y = "Geometric Mean WIS") +
#geom_point(aes(forecast_date, round(wis, digits = 2), color)), alpha = 0.05) +
facet_wrap(~ahead, nrow = 4, labeller = labeller(ahead=facets.label)) +
theme_bw() +
theme(plot.title = element_text(hjust = "center"),
legend.position = "bottom",
legend.title = element_blank()) +
scale_y_log10() +
geom_hline(yintercept = 1, size = 1.5) +
scale_color_viridis_d() +
guides(color = guide_legend(ncol = 2))
ggplotly(plot_wis, tooltip="text", height=800, width= 1000) %>%
layout(hoverlabel = list(bgcolor = "white"))
plot_wis_a <-
plot_canonical(results,
x = "ahead",
y = "wis",
aggr = GeoMean,
grp_vars = c("forecaster"),
dots = TRUE,
base_forecaster = "COVIDhub-baseline") +
labs(title = subtitle,
x = "Days ahead",
y = "Geometric Mean WIS") +
theme_bw() +
theme(plot.title = element_text(hjust = "center"),
legend.position = "bottom",
legend.title = element_blank()) +
geom_hline(yintercept = 1, size = 1.5) +
scale_y_log10() +
scale_color_viridis_d() +
guides(color = guide_legend(ncol = 2))
ggplotly(plot_wis_a, tooltip="text", height=800, width= 1000) %>%
layout(hoverlabel = list(bgcolor = "white"))
plot_wis_a <-
plot_canonical(results,
x = "ahead",
y = "wis",
aggr = Mean,
grp_vars = c("forecaster"),
dots = TRUE,
base_forecaster = "COVIDhub-baseline") +
labs(title = subtitle,
x = "Days ahead",
y = "Geometric Mean WIS") +
theme_bw() +
theme(plot.title = element_text(hjust = "center"),
legend.position = "bottom",
legend.title = element_blank()) +
geom_hline(yintercept = 1, size = 1.5) +
scale_y_log10() +
scale_color_viridis_d() +
guides(color = guide_legend(ncol = 2))
ggplotly(plot_wis_a, tooltip="text", height=800, width= 1000) %>%
layout(hoverlabel = list(bgcolor = "white"))
plot_cov80 <-
plot_canonical(results,
x = "forecast_date",
y = "cov_80",
aggr = mean,
grp_vars = c("forecaster","ahead"),
facet_rows = "ahead",
dots = FALSE) +
labs(title = subtitle, x= "Forecast date", y = "Mean Coverage 80") +
facet_wrap(~ahead, nrow = 4, labeller = labeller(ahead = facets.label)) +
theme_bw() +
theme(plot.title = element_text(hjust = "center"),
legend.position = "bottom",
legend.title = element_blank()) +
scale_color_viridis_d() +
geom_hline(yintercept = 0.8, size = 1.5) +
guides(color = guide_legend(ncol = 2))
ggplotly(plot_cov80, tooltip="text", height=800, width=1000)
plot_cov80_a <-
plot_canonical(results,
x = "ahead",
y = "cov_80",
aggr = mean,
grp_vars = "forecaster",
dots = TRUE) +
labs(title = subtitle, x= "Days ahead", y = "Mean Coverage 80") +
theme_bw() +
theme(plot.title = element_text(hjust = "center"),
legend.position = "bottom",
legend.title = element_blank()) +
scale_color_viridis_d() +
geom_hline(yintercept = 0.8, size = 1.5) +
guides(color = guide_legend(ncol = 2))
ggplotly(plot_cov80_a, tooltip="text", height=800, width=1000)
library(sf)
results_intersect <- intersect_averagers(
scores, c("forecaster"), c("forecast_date", "geo_value")) %>%
select(c("ahead", "geo_value", "forecaster","forecast_date", "data_source", "signal","target_end_date","incidence_period","actual","wis","ae","cov_80"))
kl <- function(q, p = .8) p*log(p/q) + (1-p)*log((1-p)/(1-q))
maps <- results_intersect %>%
group_by(geo_value, forecaster) %>%
summarise(wis = Mean(wis),
cov_80 = Mean(cov_80),
kl_80 = kl(cov_80)) %>%
left_join(animalia::state_population, by = "geo_value") %>%
mutate(wis = wis / population * 1e5) %>%
pivot_longer(c("wis", "cov_80", "kl_80"), names_to = "score") %>%
group_by(score) %>%
mutate(time_value = Sys.Date(),
max = max(value), min = min(value)) %>%
group_by(forecaster, .add = TRUE)
keys <- maps %>% group_keys()
maps <- maps %>% group_split()
levs <- levels(maps[[1]]$score)
# for county prediction, set geo_type = "county"
maps <- purrr::map(maps,
~as.covidcast_signal(
.x, signal = .x$score[1],
data_source = .x$forecaster[1],
geo_type = "state"))
maps <- purrr::map2(
maps, keys$score,
~plot(.x,
choro_col = scales::viridis_pal()(3),
range = switch(.y,
cov_80 = c(0,1),
wis = c(.x$min[1], .x$max[1]),
kl_80 = c(0, .x$max[1]))))
cowplot::plot_grid(plotlist = maps[keys$score == "wis"], ncol = 3)
cowplot::plot_grid(plotlist = maps[keys$score == "cov_80"], ncol = 3)
cowplot::plot_grid(plotlist = maps[keys$score == "kl_80"], ncol = 3)